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Abstract: In this paper we study the phenomenology of stars and galaxies in massive 
bigravity. We give parameter conditions for the existence of viable star solutions when the 
radius of the star is much smaller than the Compton wavelength of the graviton. If these 
parameter conditions are not met, we constrain the ratio between the coupling constants of 
the two metrics, in order to give viable conditions for e.g. neutron stars. For galaxies, we 
put constraints on both the Compton wavelength of the graviton and the conformal factor 
and coupling constants of the two metrics. The relationship between black holes and stars, 
and whether the former can be formed from the latter, is discussed. We argue that the 
different asymptotic structure of stars and black holes makes it unlikely that black holes 
form from the gravitational collapse of stars in massive bigravity. 
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1 Introduction 

The Hassan-Rosen theory (also called massive bigravity or bimetric theory) - which is the 
ghost free, non-linear theory of two interacting spin-2 fields - has garnished a lot of atten¬ 
tion concerning its phenomenological applications since it was introduced in Refs. [1, 2] in 
2012. In particular, studies of cosmological expansion histories [3-8], structure formation 
[9-23], tensor modes [24-27], integrated Sachs-Wolfe effect [20] and galactic lensing [28] 
have been performed. The possibility to couple both fields that are present in the theory 
to matter has been explored in Refs. [29-41]. In this paper, we continue the phenomeno¬ 
logical investigations of the Hassan-Rosen theory through a study of static and spherically 
symmetric (SSS) spacetimes describing, e.g., galaxies, stars and black holes. 

Since the Hassan-Rosen theory contains two metrics, the SSS spacetimes necessarily 
become more involved. In this paper we are interested in metrics that are asymptotically 
flat. These spacetimes are asymptotically classified according to the relative strength of 
the massive and massless spin-2 mode that the theory contains [42], and the conformal 
relationship between the two metrics at infinity. Concerning the black hole solutions, it 
was shown in Ref. [43] that if one assumes non-singular solutions, the two metrics must 
share a common Killing horizon. This means that black hole solutions are highly restricted. 
For star solutions, one has the option of how to couple matter to the two metrics. In this 
paper we opt for the commonly chosen approach of coupling only one of the metrics to 
matter. The theory predicts that including a gravitational source gives rise to a fixed 
relationship between the asymptotic massive and massless spin-2 modes [28, 44], In this 
paper we show that this relationship is not the same as that for black holes. This makes 
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it unlikely that the black holes that the theory contains are end-states of the gravitational 
collapse of matter. A possible cause is that the symmetry between the two metrics that 
the black holes display through the common Killing horizon, is broken when one couples 
only one of the metrics to matter. 

Spherically symmetric systems in the context of the Hassan-Rosen theory were first 
studied in Ref. [45], where, in particular, the perturbative solutions to the equations of mo¬ 
tion was published. Ref. [46] performed an extensive numerical study, and gave conditions 
for the existence of asymptotically flat black hole solutions. These solutions were further 
studied in depth in Ref. [47]. Star solutions and the so-called Vainshtein mechanism, de¬ 
scribed further below, was studied in Ref. [44], This reference is central to the analysis 
performed in this paper. Solutions for charged black holes were found in Ref. [48], and for 
rotating black holes in Ref. [49] . Stability properties of the black holes were investigated in 
Refs. [50-54], A general review of black holes in massive bigravity can be found in Ref. [55]. 

The goal of this paper is to investigate what conditions on the parameters of the theory 
that give rise to phenomenologically viable SSS solutions. Allowing for general parameter 
values by approaching the regime where the Hassan-Rosen theory becomes equivalent to 
general relativity, we constrain the ratio between the coupling constants of the two metrics. 
Furthermore, we analyse the relationship between black hole and star solutions, in order 
to see whether the gravitational collapse of stars can lead to the black hole solutions of 
massive bigravity. 

This paper is organized as follows. In Section 2 we introduce the Hassan-Rosen theory 
and the spacetime configuration under consideration. Section 3 describes the asymptotic 
solutions, and in Section 4 we state the solution for stars and their phenomenology. In 
Section 5 we constrain the phenomenology of galaxies. Section 6 discusses the relationship 
between stars and black holes, and if black holes in massive bigravity can be considered as 
end-states of the gravitational collapse of stars. We conclude in Section 7. 

We use units where G = c = 1 and M 2 = (87r) _1 . 

2 Setup 

The Lagrangian for the Hassan-Rosen theory is given by 

M 2 __ Mf , _ 

£ = -^V-det gR g -— det fR f 

4 

+ m 4 ^/- detff y^/3 n e n (Vff" 1 /) + -yZ-det g£ m , (2.1) 

n =0 

where £ m is the matter Lagrangian and e n are the elementary symmetric polynomials 
presented e.g. in Ref. [56]. Varying the Lagrangian yields the equations of motion 

3 1 
GU + m 2 J2 (“I )”= W T^, (2.2) 

n =0 J y 

+ (-l) n /3 4 —n UxYfc )v (v7^) = 0. (2.3) 

n=0 
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Here, we have defined 


2 


k = 


W n 


(2.4) 


and the matrices Y n are given in Ref. [56]. The parameter k is in principle redundant, since 
it can be put to unity through a rescaling of and the (3 n (see e.g. Refs. [56, 57]). We 
will keep it explicit, however, since it makes the limit to general relativity manifest. 

For the fields g lw and we use the following spherically symmetric and diagonal 


ansatz 


l 


ds 1 2 g = - Q 2 dt 2 + N~ 2 dr 2 + r 2 dn 2 , 

U' 2 

ds 2 = — a 2 dt 2 + -yfdr 2 + U 2 dQ 2 , (2.5) 

where a prime signifies a derivative with respect to r. This form for g gv and f gu is the 
most general diagonal form of the metrics after using the possibility of doing a rescaling of 
the radial coordinate. Notice that f /w can equivalently be written 

ds 2 = - a 2 dt 2 + Y~ 2 dU 2 + U 2 dn 2 , (2.6) 

and U(r ) be interpreted as the radial coordinate for the /-metric. 

The energy density and pressure are given by p(r ) = — Tq and P(r) = T -/3 (summation 
over i implied), and they satisfy the following conservation equation: 

P' = ~^(P + p). (2.7) 

In this paper, we will combine analytic and numerical studies. For the numerical 
analysis, we follow Ref. [46] and put the equations of motion in the following form: 


N' 

= Pi (r 

, Q, N, Y, 

U, p, P, c. 

,m 2 ,^i,/32,/3 3 

,k) > 

Y' = 

= Pi (r, 

Q,n,y, 

U, p, P, c, 

m 2 ,/3i,/3 2 ,/33, 

k) , 

U' ~- 

= P'i (r, 

Q,n,y , 

u , P, P, C, 

m 2 ,/3i,/3 2 ,/3 3 , 

«) , 

Q'~- 

= P 1 {r, 

Q,n,y, 

U, p, P, c, 

m 2 1 Pi, P 2 , P 3 , 

k) , 

p’ -- 

= P 5 {r, 

Q,n,y, 

U, p, P, c, 

"i 2 ,^i,/3 2 ,/3 3 , 

k) , 


where c is defined below. The function a can be solved for directly once the other fields are 
given. When p = P = 0, i.e. in vacuum, J -3 become independent of Q. In vacuum, 

one thus first solves three first order equations for N, Y and U, and then integrate J -4 to 
get Q. When p and P are non-vanishing, the five first order differential equations instead 
have to be solved simultaneously. 

1 For a non-diagonal ansatz, the equations of motion constrain the solution to be identical to the 

Schwarzschild-AdS/dS metric, as shown in Ref. [45, 46]. 
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3 Asymptotic structure 

Since we are interested in solutions that are asymptotically flat, the metrics should approach 


9/iv ^ 


I [IV 




(3.1) 


at infinity. Here c is an asymptotic conformal factor between the two metrics. In order for 
Eq. 3.1 to be an solution, we need to impose 


Po =-3/3ic-3/? 2C 2 -/3 3 c 3 , (3.2) 

Pi = ~ Pic ~ 3 - 3/3 2 cE 2 - 3/3 3 c“ 1 , (3.3) 


to cancel the cosmological constant terms for g^ and f^ 2 

Linearizing around the flat space backgrounds, i.e. expanding the metric components 
as Q = 1 + 5Q, N = 1 + 5N, a = c(l + 5a), U = cr (1 + 5U ) and Y = 1 + 5Y, gives 


5Q 

5N 

5a 

5Y 

5U 


C\ ( C 2 kc 2 _ r 
2 r 


= - — + e -m 9 r 


_Ci C 2 kc 2 (1 + m g r) r 
2r 2 r 


9l _ 

2 r r 

Cl _ C 2 (1 + m g r) ~-m n r 
2 r 2 r 


(l + kc 2 ) C 2 (l + m g r + m 2 r 2 ) _ mgT 
2 m^r 3 


(3.4) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 


These solutions, first appearing in Ref. [45], are well-known and have been presented on 
several occasions in the literature. The parameters Ci and C 2 regulate the strength of the 
massive and massless modes. The graviton mass m g is given by 


2 2 
m g = m 


1 + 


KC* 


(Pic + 2 p 2 c 2 + p 3 c 3 ) 


(3.9) 


Let us discuss the free parameters that we have at our disposal. Of the five p n , two 
have been fixed in order to yield asymptotically flat solutions. This leaves Pi, p 2 and Ps as 
free theory parameters, m 2 is not a free parameter since it can be absorbed into the p n : s. 
We will keep it explicit, however, since it sets an overall length scale when the p n :s are of 
order unity. As mentioned above, k is also redundant since it can be put to unity through 
a rescaling of and the p n :s. Since it is important for discerning solutions that lie close 
to those of general relativity, we will, however, keep it explicit. Added to this, we have the 
conformal factor c. On the whole, for vacuum solutions, we have four global parameters, 

2 Notice that in Refs. [46] and [47] the parametrization 

/3n = (-l) n+1 Q (3 - n) (4 - n) - (4 - n) a 3 - a^j 

was used, for which c = 1. 


- 4 - 







the three /3,:s and c, together with the local parameters C\ and C 2 , which controls the 
strength of the massless and massive modes. As discussed later, including a gravitational 
source fixes the relation between C 1 and C 2 . 

The equation of motion have the property that under the rescaling 


N(r) —> N(Xr), Y (r)—*■ Y (Xr), Q(r) —>• Q(Xr), a{r) —> a(Xr), 

f) ~P TY1 

U(r) —U(Xr), pP (3.10) 

a solution is mapped onto a new solution [46]. We will interchangeable use ry (defined 
below) or X g = m” 1 as radial coordinate. 

The linear solutions are valid up to the radius where higher order terms become im¬ 
portant. This radius is usually called the Vainshtein radius, and was first identified in 
Ref. [58] in 1972. In massive bigravity, the Vainshtein radius is 


tv = 



(3.11) 


where M to t is defined as the total mass of a source. In Ref. [58], Vainshtein also conjec¬ 
tured that there should exist a mechanism, later dubbed the Vainshtein mechanism, that 
effectively restores general relativity inside the Vainshtein radius. That this exists in the 
context of massive bigravity for SSS spacetimes was shown in Ref. [46] for the case of 
k —} 0, and in Ref. [44] for the r <C X g limit. It is important to note, however, that the 
existence of the Vainshtein mechanism depends on the specific choice of the parameters. 

For recent phenomenology concerning the Vainshtein mechanism, see Refs. [23, 59, 60] 
and references therein. 


4 Stars 


In this section we study the phenomenology of stars in massive bigravity. As a source, we 
use a star with constant energy density p*, pressure P(r) and radius r*. The pressure has 
to satisfy the conservation equation (2.7), and vanish at the surface of the star. The mass 
interior to r is 

M ( r ) = J o P(r)r 2 (lr, (4.1) 

and the total mass of the star is thus 


M _ £*[*_ 

* 6M„ 2 ' 


(4.2) 


We have three effective scales for the stars: r*, ry and X g . We will assume that r* « X g , 
and comment on both the r* < ry scenario as well as r* > ry. 

As shown in Refs. [28, 44], the introduction of a source fixes the relation of C\ and C 2 
in the linear solutions to 


c 2 = - 


2C\ 


(4.3) 
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and C i = 2Mtot/(l + kc 2 ). The linear solutions then become 


5Q = 
5N = 
5a = 
5Y = 
5U = 


Mi 


tot 


4kc 2 _ 
1 + —-—e mgr 


1 + 


3 

2 kc 2 


r(l + kc 2 ) 

Mtot 

r( 1 + kc 2 ) 

Mtot 

r( 1 + kc 2 ) 

M t °t 2 _ 

- 77 —-oT 1 — o (1 + m g r)e 

r{ 1 + kc 2 ) [3 

2M tot [1 + m 9 r + {m g r) 2 ]e~ rnar 

3 r{m g r) 2 


1 _ 2p~ m 9 r 

3 


(4.4) 

(4.5) 

(4.6) 

(4.7) 

(4.8) 


Asymptotically, the fields thus look like a massless general relativity (GR) like term plus 
a Yukawa term. They exhibit the usual vDVZ-discontinuity [61-63] which can be probed 
observationally. As r 2 > X g , the Yukawa term decays, however, and the fields look identical 
to general relativity. Is is only when r <a 9 , or when higher order terms become important, 
that we can expect any observational signatures. 

When massive bigravity is used for cosmological applications, for kc 2 ~ 1, we expect 
A g to be of the order of the Hubble scale today. It is then an excellent approximation 
that r* <C X g . For this framework, it was shown in Ref. [44] that it is possible to obtain 
approximative analytical solutions by assuming that all fields and their derivatives are close 
to the flat space background, with the exception of U/r. Defining 3 


_ U 

V =-1, 

cr 


a = — 


P = 


foe 2 + /V 
foe + 2foe 2 + foe 3 ’ 
foe 3 


foe+ 2 foe 2 + foe 3 ’ 
the metric perturbations can be expressed as 


rW P ( r ) r 2 M (r) KC 2 m 2 r 2 

rSQ' = fo i n + 9 


5N = — 
rSa' = 


2Mg 

M(r) 


2(1 + KC 2 ) 


A* 


P 


;fo 


2 2 2 
nc m g r 


r 2(1 + kc 2 ) 
m 2 r 2 (r + rfo' 


5Y = 


2(1 + kc 2 ) (1 + nY 

2 2 
m g r 

2 (1 + /i) (1 + KC 2 ) 


2 P 3 

H - afo + -fo 
1 


/i + 2 fo 2 + - (2 - 2a - fo /r 3 

O 


P + (1 — a ) + g (1 — a + fo) M 3 


(4.9) 

(4.10) 

(4.11) 

(4.12) 

(4.13) 


3 Note that the definitions of a and f) are generalized compared to Ref. [44] in which /3ic + 2/32c 2 + [5z(? 

was normalized to unity. 
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K = C = 1 K ~ C = 1 




Figure 1 . Le/£ panel: Allowed region (hatched) for a and /3 using k = c = 1. Right panel: Allowed 
region (hatched) for fa and fa using a normalization where fa c + 2 fac 2 + fac 3 = 1 and n = c = 1 . 


These fields are thus functions of M(r), P(r) and //, where // satisfies a seventh-degree 
polynomial: 


3(l + kc 2 ) /t + 6 (l + kc 2 ) (1 — a) // 2 + 
i [6(l + kc 2 ) a 2 — 2 (17 + 18kc 2 ) cc + 4 (l + kc 2 ) 0 + 10 + 9kc 2 ] /z 3 + 

O 

^ [6 (l + kc 2 ) a 2 — (7 + 9kc 2 ) a + 4 (l + kc 2 ) 0 + l] /z 4 + 
^[2(l + 3kc 2 ) a 2 — (l + kc 2 ) (3 2 + 2 (l + 2 kc 2 ) 0 — 4 a0 — 2a] // 5 

O 


2 2o2 6 1 2o2 7 

--KC 0 /i - -kc 0 n 


(l + KC 2 ) (1 + fj) 


2 r 


m- 


2M (r) 


(1 - 0 // 2 ) - (1 - 2 a// + 0 /z 2 ) 

9 


(4.14) 


The function // satisfies — 1 < // < 0 for all physically relevant cases. 

In Sec. A, we show that real valued solutions to Eq. 4.14 that approach zero at infinity 
(which corresponds to the asymptotically flat solutions) exist if a > —\/y/~j3. Furthermore, 
one must also have a < —d\/d 2 when d 2 < 0 , where 


d\ = 1 + 3 kc 2 — 6-\/0(l + kc 2 ) + 3/3(1 + kc 2 ), 

d 2 = —1 + 6y//3(l + kc 2 )(1 + 0) - 0(13 + 12kc 2 ). (4.15) 


These constraints are depicted in Fig. 1 and are more restrictive than those presented in 
Ref. [44], In terms of fa, we can write a > and 0 > 1 as (using the normalization 

where fac + 2fac 2 + fac 3 = 1 ) 

fac 3 > 1 , 

fac 2 < y/ fac 3 - fac 3 < 0 , 

fac = 1 - 2fac 2 - fac 3 > (y/ fac 3 - l ) 2 > 0. (4-16) 
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That is, we need /3 2 to be strictly negative and /3i and /I 3 to be strictly positive. 

For the phenomenological analysis, we will use the following definitions of the poten¬ 
tials: 


4> = Q 2 - 1 


2(5Q, 


^ = 1 - N~ 2 ~ 2(5 Ah 


From Eqs. 4.10, 4.11 and 4.17, we have 

2 M (r) 


T = -- 


2 2 2 
m g ncr 

1 + KC 2 


2 

— oy -|—— 


2 P(r) 2 2 M(r) 

v =r * = -Mi r -^—r 1 


2 2 2 
rrigKcr 

1 + KC 2 




/V 


(4.17) 

(4.18) 

(4.19) 


where v is the circular velocity. From this we see that as long as // stays real and finite as 
r —> 0, we will recover GR at small radii, as long as the potentials are small. 

Outside the source, we have 

3 

K.rr / r \ 

1 + 


T = -- 


2 Mi 


tot 


r^>ry 


r 

2 M tot 


r<F = 


>- 

r 

2 M to t 


1 - 


1 + KC 2 
KC 2 

3 (1 + KC 2 ) 


r 

ry 


2 , PP 
[i — an -|—— 


(4.20) 


r~^>ry 2M to t 


1 - 


1 + 


KC 


1 M 

1 + kc z \ryJ 


KC 

3(1 + KC 2 ) 






(4.21) 


(4.22) 


The limit r ry is derived by noting that we can neglect higher order terms in /U, and 
solve Eq. 4.14 as 

1 < TV ^ (4.23) 


V = ~o l - / 

3 V r / 

Note that this “asymptotic” value is only valid far inside the Compton wavelength of 
the graviton and represents the maximal deviation we expect from GR. The deviation 
is monotonically increasing with kc 2 and has a maximal value of 1/3. For r > A 9 , we 
recover GR again. If we are far outside the graviton Compton wavelength, the exponential 
term will be negligible and GR is recovered. If we are far inside the Vainshtein radius 
ry ~ (MtotA 2 ) 1 / 3 , and GR is restored again. We thus expect the largest deviations from 
GR to happen at ry < r < \ g . As an example, assuming m g = Hq, for the Sun, we have 


r-O 

V 


4 • 10 18 m « 140 pc rs 2.7 • 10 7 AU. 


(4.24) 


We thus only expect the gravitational field from the Sun to be modified on scales much 
larger than the distances to its closest star neighbours, where they of course are completely 
negligible anyway. In the solar system (r ~ 1 AU), deviations are of order (r/ry ) 3 ~ 10~ 21 . 
This value is way below current observational constraints showing that on AU scales, 
deviations from the inverse square force law is < 10 ~ 9 [64], This observational constraint 
indicate that X g > 1 kpc, except for the case of kc 2 <C 1 when constraints on X g will be 
weaker. 
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For parameter values not fulfilling the requirements given, we may still have everywhere 
real solutions if the radius of the source is bigger than its Vainshtein radius. For example, if 
the source has a constant density (and zero pressure), p will be constant inside the source. 
As shown in Ref. [23], the requirement of having sources larger than their Vainshtein 
radii corresponds to them having densities smaller than order of the critical density of the 
Universe, if m g ~ m ~ Hq. We still might be able to have more compact sources if we let 
k be very small since m g oc m,K~ 1 ^ 2 . For a source with mean density p*, for small k, we 
can write Eq. 4.14 at the surface of the source 

3p + 6 (1 — a) p 2 + \ [6a 2 - 34a + 4/3 + 10] p 3 4 - 

O 

- [ 6 a 2 — 7a 4- 4/3 4-1] p 4 4— [ 2 a 2 — /3 2 4- 2/3 — 4a/3 — 2a] p 5 
3 3 

(4 ' 25) 

where the critical background density of the Universe today is given by 

p CI = 3M 2 Hq ~ 1.88 • 1(T 29 h 2 g cm~ 3 . (4.26) 

Demanding that solutions exist down to the surface of a neutron star for which 

Pneutr ~ 4 • 10 14 gcm” 3 , (4.27) 

we generally need (p*/p cr ) {Ho/m g ) 2 « (p*/p cr ) k (Ho/m) 2 to be smaller than order one, 
which means that k should be less than the ratio of the critical density of the universe and 
the source density, assuming m ~ i7o - 4 For general values of m we get 

« < 10 ” 44 . (4.28) 

Alternatively, we can constrain the Compton wavelength of the graviton 

< \l ——— r H — 28km. (4.29) 

V Pneutr 

where rn = H^ 1 ~ 1.3 • 10 26 m. Note however that this very restrictive limit only needs 
to be fulfilled for parameter values not fulfilling the ones illustrated in Fig. 1. 

In Ref. [65] a limit of y/R < 10 -1 ' was derived in order to push scalar instabilities 
back before BBN. For this to work, we need at least two /3* 7 ^ 0. For background and 
perturbation solutions, the main idea is that in the limit of At —>• 0, the ratio between the 
scale factors of the two metrics goes to a constant determined by the values of the /3 i and 
c. This gives a cosmological constant-like contribution to the Friedmann equation and 
well-behaved perturbation theory. 

4 Although fields are not weak at the surface of a neutron star, using the polynomial equation is sufficient 

to obtain order-of-magnitude estmates. 
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Figure 2. Deviations from the general relativity predictions for the potential T and the circular 
velocity v 2 for the case of n = c = l = 1 (where l = r g /rv ), a = 1 and /? = 4. 


5 Galaxies 


We now turn to the phenomenology of galaxies. In the dark matter paradigm, there is 
a somewhat unexpected large correlation between the distribution of baryonic and dark 
matter. One of the main arguments for MOND is that it is able to explain this correlation 
on galactic scales [66]. However, it fails on larger scales [67]. The Vainshtein radius on the 
other hand naturally adapts to the scale of the object. 

We add a galactic source (with negligible pressure) with density profile 

p{r) = p 0 r~ q , (5.1) 


being truncated at r = r g = lry 1 where the parameter l sets the compactness of the galaxy 
and is of order one, or slightly lower, for m g ~ Hq. We then have 


M(r) =p 0 


Mtot — po 


q 

2(3-g)’ 

3 —q 
r 9 

2(3-g)' 


(5.2) 


We can now write (inside the galaxy; outside the galaxy the solution is given by Eqs. 4.20) 


2 M(r) 


1 + l 3 ~ q 


KC 


1 + KC 2 


(£) " (" - 0:1,2 + 



V 


2 


r4>' = 


2 M(r) 

r 



KC 


1 + KC 2 



(5.3) 


This typically gives result like in Fig. 2 where k = c = l = 1, a = 1 and (3 = 4. 

The observed gravitational lensing and dynamical properties of elliptical galaxies are 
consistent with general relativity predictions, to an accuracy of ~ 5% [28, 68]. We have 
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three ways to make our model consistent with lensing constraints. The first is to make the 
Compton wavelength so small that we are well outside it for the lensing and dynamical 
observations (basically, the velocity dispersion of stars). The lensing radii typically are 
~ 5 kpc and the velocity dispersion integrated out to similiar radii. In order not to be in 
conflict with the observed constraints, we thus need X g < 0.5 kpc. However, as noted in 
Sec. 4, such small values are ruled out by Solar system constraints. 

The second possibility is that the so called gravitational slip 7 - the ratio of the 
gravitational potentials experienced by massive and massless particles - is small. The 
largest deviations from general relativity predictions are found between the Vainshtein 
radius and the Compton wavelength where 

7 = __ (54) 

7 3(1 + 3kc 2 ) 1 ’ 

Using data from the strong gravitational lens sample observed with the Hubble Space 
Telescope Advanced Camera for Surveys by the Sloan Lens ACS (SLACS) Survey [69], we 
constrain kc 2 < 0.1 at 2 a. 

The third possibility, valid for parameter values for which we have a functioning Vain¬ 
shtein mechanism, is to make sure that we are well inside the Vainshtein radius, 


ry 



M 


10 11 M 0 


1/3 



kpc. 


(5.5) 


For large kc 2 , we typically need to be a factor of 10 inside the Vainshtein radius not to be 
in conflict with observational limits, corresponding to m g /H q < 40 or A 5 > 0.1 Gpc. 

To summarize, strong lensing galaxy systems constrain the graviton Compton wave¬ 
length X g to be either smaller than ~ 0.5 kpc or larger than ~ 0.1 Gpc, or the combination 
kc 2 to be smaller than ~ 0.1. However, X g < 0.5 kpc is disfavoured by Solar system con¬ 
straints. We also note that we generally expect the velocity dispersion in galaxies and 
galaxy clusters to increase as compared to the general relativity prediction, on scales sim¬ 
ilar to the sizes of the systems if m g ~ Hq or slightly larger. This will have an effect 
on the predicted abundance of dark matter in these systems, namely that we need less 
dark matter than in the general relativity case. However, since we maximally expect the 
velocity dispersion squared to increase by a factor of 1/3, the effect is not large enough to 
completely evade the need for dark matter in galaxies and galaxy clusters. 


6 Vacuum solutions 

In the previous section, we studied stars, galaxies and their phenomenology. In this section 
we comment on the relationship between the star solutions and vacuum solutions, such as 
black holes. Our chief interest here is to understand if the bimetric black holes can be the 
end-state of the gravitational collapse of massive stars. 

Vacuum solutions in massive bigravity were studied extensively in Ref. [46] . Following 
the proof of Ref. [43]— that for non-singular metrics there has to be a common Killing 
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Figure 3. Left panel: The field N for different choices of c, for a black hole with rn/Xg = 0.04 and 
/?2 — /3s — 0. The fields approach the Schwarzschild solution close to rjj- Varying the parameter c 
shows that several possible solutions exist for a given rn/X. Right panel: The constant u, given by 
JJ/r as r —> r h, as a function of c, showing a close to linear relationship. 

horizon—we expand the fields N, Y and U close to the horizon, situated at r = rh, as 
N 2 = ^2 a n (r - r h ) n , Y 2 = ^ b n (r — r h ) n , U = ur h + ^ c n (r - r h ) n . (6.1) 

n>1 n >1 n >1 

From Eqs. 2.8, the coefficients a n , b n and c n can all be expressed in terms of u and a±, 
where u is arbitrary and a\ satisfies a quadratic polynomial with coefficients depending 
on u and the parameters of the theory (i.e. c and the /% parameters). Since there are 
three equations of motion, and three free parameters (u, C\ and C 2 ) there exists at most 
a discrete set of solutions for a given value of c and the /3j parameters. The structure of 
these solutions was investigated extensively in Ref. [47], in the case of c = 1. This was 
done through a shooting method, where u, C\ and C 2 were varied until the solution with 
asymptotic flatness was found. In this paper we have performed a similar numerical study, 
but with general c. Our results are in agreement with Refs. [47] and [46] wherever they 
overlap. 

It was found in Ref. [47], that for a given value of the /% parameters, the solutions are 
classified by r^/Xg, i.e the ratio between the horizon and the Compton wavelength of the 
graviton. An upper bound for rh/X g is 0.876, a value related to the Gregory-Laflamme 
instability (see Ref. [70] for an interesting discussion of this result). Above that bound, 
only the bi-Schwarzschild solution exists (i.e. g is equal to the Schwarzschild solution, 
and fpu = (?g t ,, v )- The minimum value of rhjX g depends on the model under consideration. 
The conjectured parameter structure presented in Ref. [47] is that when ^3 is non-zero, 
solutions cease to exist below a critical value of Vh/X g (which excludes realistic astrophysical 
black holes). When /S 3 = 0, /?2 > 1 and f3\ < —1, black hole solutions exist for all values of 
rh/Xg below the Gregory-Laflamme bound. 

Moving beyond the case of c = 1, we show in Fig. 3 the metric field N for different 
values of c and with ru/X g = 0.04. This shows that several possible black hole solutions 
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Figure 4. Left panel: The function U/r solved using the full equations of motion numerically 
(dashed) and using the approximate solution, given by Eq. 4.14 (dotted). U/r departs from the 
constant solution predicted by the approximate solution when the other metric fields become non¬ 
linear. Right panel: The metric function N divided by the GR-solution and N divided by Q. For 
the GR-solution, the Schwarzschild radius is given by rs/X g = 10 -4 ; this ensures that the horizon 
of the GR and bigravity solutions conicide. In both the left and right panel, C\/\ g = 5 x 10 -5 , 
C 2 = —2/3 x Ci and (3\ = 7, /3 2 = —5, (3s = 4, c = n = 1 (these specific values ensure that the 
solution exists within the Vainshtein radius). N/Nqr and N/Q approach unity as C\/\ g decreases. 


are possible for fixed rjy/A 9 , as long as c is varied. We also display the relationship between 
the constants u and c. 

Stars and black holes. Concerning the relationship between the star and black hole 
solutions, we note the following: First of all, u < c for the stars, but u > c for the black 
holes. Secondly, for star solutions to exist inside the Vainshtein radius, we must have 
/? 3 C 3 > 1. For the black holes, we must instead have /?3 = 0 for solutions to exist for all 
rh/^g < 0.876, according to the conjecture of Ref. [47]. Finally, the asymptotic structure is 
different as compared to the black holes and stars. For stars, we have C 2 /C 1 = —2/3. For 
the black holes, while a full parameter scan is beyond the scope of this paper, we conjecture 
that all black holes satisfy 


0 < 


C 2 c 2 

Ci 



( 6 . 2 ) 


This conjecture follows from a numerical analysis, where we find that the point C^c 2 = 
2C\/?> (in the following we put k = 1) marks a transition for the behaviour of N. Above 
this value, i.e. C 2 c 2 > 2Ci/3, N will generically become larger than unity. 5 Below this 
value, N will become less than unity. Furthermore, for C 2 > 0, U/r will grow larger than c 
as one integrates from infinity towards lower r, and for C 2 < 0, it will become smaller. The 
point C 2 = 0 corresponds to the Schwarzschild solution, and as C 2 —X 0, r^/ \ g approaches 
the value 0.876 given by the Gregory-Laflamme instability. For the black holes, we have 


5 As a sideremark, we note that for C 2 C 2 > 2Ci/3, there exist solutions where all the metric fields beside 
a go like ~ 1/%/r, inside the Vainshtein radius, for the /3 2 = /?3 = 0 model. The implication of these 
solutions will be investigated in an upcoming work. 
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that N should become less than unity (and eventually approach zero), and U/r should 
be larger than c. Thus, for the black holes, we should have 0 < C 2 C 2 < 2Ci/3. This 
is also confirmed for the black hole solutions that we have studied. There thus seem to 
be qualitative difference concerning the overall sign of the massive Yukawa modes when 
comparing stars and black holes. This stands in contrast to the case of general relativity, 
where a spherical collapse of a massive star into a black hole does not change the asymptotic 
spacetime structure. 6 

What happens, then, in vacuum when the asymptotic structure of stars is imposed? 
Solving the full numerical system, we find that the fields of g^ v approach the Schwarzschild 
solution (for parameters that satisfy the bounds given in Eq. 4.16). The function U/r 
remains constant in a region inside the Vainshtein radius, but starts to grow close to the 
horizon of g^ u . The fields a and Y remain small and finite. We depict this scenario in 
Fig. 4. 

An interesting curvature invariant, introduced in Ref. [43], is 

~ 2 U' 2 N 2 2 U 2 


i = frg* 1 '= qs + 


Y 2 


+ 


(6.3) 


This function remains finite for all non-singular metrics, in particular for the black hole 
and star solutions. For the vacuum solution shown in Fig. 4, it does, however, diverge 
close to the horizon of the g metric. This is related to the fact that there is no common 
horizon for both g^ u and when C 2 /C 1 = —2/3 in vacuum. 

Instabilities. Let us also discuss the instabilities that are present for the bi-Schwarzschild 
solutions. It was shown in Refs. [50, 51] that there exists unstable modes when the horizon 
radius of the source is less than the Compton wavelength of the graviton. This instability 
is, however, rather mild, with a timescale equal to the inverse graviton mass. When the 
latter is of the same size as the Hubble scale today, this means that the instability will 
require the entire lifetime of the universe to grow significantly. It does, therefore, not have 
to be important for astrophysical black holes. Intriguingly, the instability was shown to 
be absent for the non-diagonal bi-Schwarzschild solutions [53], as well as for the partially 
massless case [52]. Now, as was argued in Refs. [50, 51, 70], the instability shows that the 
bi-Schwarzschild solution can not be considered the end-state of a gravitational collapse. 
It is unclear whether the other black hole solutions, with massive hair, are stable or not. 
On the whole, then, there are two reasons why the end-state of gravitational collapse is 
unclear: the instabilities present for the bi-Schwarzschild case (which could also be present 
for the other black hole solutions), and the different asymptotic structure of stars and black 
holes. 

To summarize, there is a qualitative difference between the star and black hole solu¬ 
tions. The end state of a collapse of a star is therefore uncertain. It could lead to a novel 

6 It is possible that during the collapse process, information concerning the change of the asymptotic 
structure could propagate outwards at a finite speed. This information would take an infinite time to change 
the asymptotic structure, a process which has been observed in general relativity during the collapse of 
massive scalar fields [71]. This could potentially reconcile the different asymptotic structure for stars and 
black holes. It is an open question whether this is the case in massive bigravity (we thank the referee for 
pointing this out). 
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spherically symmetric solution that as of yet has not been discovered. It might lead to a 
time-dependent solution that does not settle down into a static final state. It seems un¬ 
likely, however, that it will lead to the black hole solution that share a common horizon for 
g and /. We therefore conjecture that black holes in massive bigravity can not be formed 
from the collapse of stars. This is probably due to the fact that the black hole solutions 
share a symmetry between g jJ:l/ and f fiu (i.e. a common horizon), whereas the coupling of 
matter to only one metric, e.g. g, breaks this symmetry. An interesting question is whether 
this conjecture also holds true when coupling matter to both fields. 

7 Conclusions 

In this paper, we have investigated the phenomenology of stars and galaxies in massive 
bigravity. Furthermore, we have discussed the relationship between black holes in massive 
bigravity and stars. 

For the stars, we have been interested in the existence of solutions where the radius 
of the star is much smaller then the Compton wavelength of the graviton. The latter is 
usually assumed to be of the order of the Hubble scale of the universe today, when massive 
bigravity is used for cosmological applications. The parameter constraint that we found, 
which generalizes earlier work in Ref. [44], states that @2 needs to be strictly negative and 
f3\ and @3 needs to be strictly positive. If these conditions are not met, we have shown 
that the ratio between the Planck masses of the two metrics needs to be less than 10“ 22 , 
when the length scale of the theory is of the order the Hubble scale today. 

Moving on to galaxies, we show that the graviton Compton wavelength \ g either has to 
be so small (less than ~0.5 kpc) so that the massive Yukawa mode does not produce sizable 
deviations between the lensing and dynamical observations. This is, however, in conflict 
with Solar system measurements. Another possibility is that \ g is so large that the galaxies 
fall within the Vainsthein radius. This requires X g >0.1 Gpc. Yet another possibility is 
that kc 2 <0.1, which makes the deviation in the gravitational slip undetectable. 

Finally, an open and interesting question, that deserves further studies, is the end-state 
of gravitational collapse. In general relativity, the asymptotic structure is unchanged as 
a star undergoes spherical collapse to a black hole (a fact related to Birkhoff’s theorem). 
In massive bigravity, we find that the asymptotic structure of stars and black holes is 
qualitatively different. This is related to the sign of the massive Yukawa mode. This 
makes it unlikely that the black hole solutions are end-states of gravitational collapse. 
This could potentially be related to the fact that for the black holes g gv and f gv have a 
common Killing horizon. This symmetry is, however, broken by stars, since only one of the 
metrics couple to matter. It would therefore be interesting to investigate the star solutions 
when coupling both metrics to matter. 

A Real solutions 

In this Appendix we derive constraints on parameter values needed to have static, spheri¬ 
cally symmetric solutions that are asymptotically flat and valid att all radii. We will make 
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numerous references to the left hand side (LHS) and right hand side (RHS) of the polyno¬ 
mial equation (4.14). Assuming that we are outside the source, M(r) = M tot , the pressure 
is zero and we start by noting that the the RHS is zero at /j = — 1 and /r = ±l/\//3 (for 
(3 > 0; for /3 < 0, the only root is at /j = — 1 ) and is being divided by r 3 . As r -> oo, 
the RHS becomes flat and as r —>• 0, RHS —>• Too, except at the points where it is zero. 
Defining (for the pressureless case) 

h = {~ ) 3 = (! + (! - /V) (A- 1 ) 

1 + KC Z \ry J v ' 


we have 


which is zero at /r = — 1 . 




= 2 (1 + //) (l — 2(3 fi — 3/5/x 2 ) , 


dfi 

Furthermore, 
d 2 h 

dji 


2 = 2 (l — 2/3 — 10/3// — 9/5/u 2 ) , 


(A.2) 


(A.3) 


which at // = — 1 is 2(1 — /?). This is negative for (3 > 1 and vice versa. 

The LHS on the other hand has a shape that is fixed by the values of /3j, //, n and c. It 
is always zero at // = 0 and 


= 3(l + kc 2 ) > 0. (A.4) 

fi =0 

Since we want our solutions to be asymptotically flat we need /r —>• 0 as r —> 00 . The 
limiting value for fj, is either /r = — 1 if (3 < 1 or fi = —l/\f]3 if [3 > 1. We first show that 
(3 < 1 is not an option. For fi = — 1, the LHS becomes 

LHS = ^(l + 2a + /3) 2 . (A.5) 

o 

This has to be smaller than or equal to zero in order for solutions not to become imaginary 
as r —>• 0, since the RHS gets arbitrarily negative close to fi = —1. This means that we 
need to set /3 = —1 — 2a. Close to fi = — 1, we can expand the RHS and LHS sides as 

LHS ~ —(1 + a) [3 + 2kc 2 + a( 3 + acc 2 )] (/j, + l) 2 , 

3 

RHS - -2(1 + kc 2 )(1 + a) ) 3 (/j + l) 2 , (A.6) 


dLHS 

dfj, 


showing that we will not have real solutions as r —>• 0 since the RHS always will be less 
than the LHS for some finite r. 

For (3 > 1, the question is whether we have a real solution for which n = [— l/y/]3, 0] 
for r = [0,oo], For /r = — l/yf]3, we can write the LHS as 


d\ = 1 + 3 kc 2 — 6 y^d(l + kc 2 ) + 3/3(1 + nc 2 ), 

d 2 = -1 + 6y^/3(l + kc 2 )(1 + /3) - /3(13 + 12kc 2 ). (A.7) 
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Figure 5. Graphic solution for k = c = 1, /3 = 4 and a = — 1 (left panel) and a = —3 (right 


panel). If the LHS has a maximum or inflection point over the interval /x = [— l/y/Ji, 0], solutions 
for which fi —> 0 as r —> oo will become complex at some finite value of r as illustrated in the right 
panel. This rules out regions where very low values of a. 

To have LHS < 0, we need 



(A.8) 


We also need the LHS not to have a maximum or inflection point over the interval /x = 
[— 0]. This rules out regions where a < —d\/d 2 when cfo > 0 as illustrated in Fig. 5. 

We are thus left with a > —\/y/jd, with the additional constraint a < —d\/d 2 for regions 
where c ?2 < 0, see Fig. 1. 
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